function [u,flag,normres,Q] = gmres(varargin)
% GMRES  Iterative solution of chebfun operator equations.
%
% U = GMRES(A,F) attempts to solve the operator equation A(U)=F, where F 
% and U are chebfuns and A is a function handle defining a linear operator 
% on chebfuns.
%
% U = GMRES(A,F,RESTART) chooses a restart parameter. Use Inf or [] for no
% restarts. The default is Inf.
%
% U = GMRES(A,F,RESTART,TOL) tries to find a solution for which the
% residual norm is less than TOL relative to the norm of F. The default
% value is 1e-10.
%
% U = GMRES(A,F,RESTART,TOL,MAXIT) stops after MAXIT total iterations. This
% defaults to 36.
%
% U = GMRES(A,F,RESTART,TOL,MAXIT,M1INV,M2INV) applies the left and right
% preconditioners M1INV and M2INV, defined as functions. They should
% approximate the inverse of A. Note that this usage of preconditioners
% differs from that in the built-in GMRES.
%
% [U,FLAG] = GMRES(A,F,...) also returns a convergence FLAG:
%     0 GMRES converged to the desired tolerance TOL within MAXIT iterations.
%     1 GMRES iterated MAXIT times but did not converge.
%
% [U,FLAG,NORMRES] = GMRES(A,F,...) also returns a vector of the relative
% residual norms for all iterations. Note the output ordering is not the 
% same as for built-in GMRES. This calling sequence will also print out
% updates on the progress of the iteration.
%
% See also GMRES, CHEBOP/MLDIVIDE.
%
% EXAMPLE
%
%   % To solve a simple Volterra integral equation:
%   f = chebfun('exp(-4*x.^2)',[-1 1]);
%   A = @(u) cumsum(u) + 20*u;
%   u = gmres(A,f,Inf,1e-14);

% Copyright 2011 by The University of Oxford and The Chebfun Developers. 
% See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.

% Parse inputs and supply defaults.
defaults = {[],[],Inf,1e-10,36,[],[]};
idx = nargin+1:length(defaults);
args = [varargin defaults(idx)];
[L,f,m,tol,maxiter,M1inv,M2inv] = deal( args{:} );
showtrace = (nargout>2);

if m==0, m = Inf; end;                           % no restarts
m = min(m,maxiter+1);  % avoid warning about Inf in for loop

u = chebfun(0,domain(f));
normb = norm(f);
r = f;

if ~isempty(M1inv), r=M1inv(r); end
normres(1) = norm(r)/normb;
j = 1;                                           % total iterations
while (normres(j) > tol) && (j<maxiter)          % outer iterations
  H = [];
  QTb = norm(r);
  Q = r/QTb ;                                    % Krylov basis
  for n = 1:m                                    % inner iterations
    % Next Krylov vector, with preconditioners.
    q = Q(:,n);                                    
    if ~isempty(M2inv), q = M2inv(q); end       
    v = L(q);
    % KLUDGE: The auto-diff facility will overwhelm memory if not shut down.
    if ~isempty(M1inv), v = M1inv(v); end       
    % Modified Gram-Schmidt iteration.
    for k = 1:n                     
      H(k,n) = Q(:,k)' * v ;             
      v = v - H(k,n)*Q(:,k);  
      v.jacobian = anon([],[],[],Inf);  % disable auto-diff
    end 
    H(n+1,n) = norm(v);                          
    Q(:,n+1) = v / H(n+1,n);                     % new basis vector
    
    % Use QR factorization to find the residual norm.
    % TODO: This could be made more efficient (worthwhile?).
    QTb(n+1,1) = 0;                              % by orthogonality
    j = j + 1;
    [P,R] = qr(H);  
    normres(j) = abs(P(1,n+1)*QTb(1))/normb;
    
    % Done?
    if normres(j) < tol
      showtrace = false;
      flag = 0;
      break
    end
    
    % Trace iteration--for debugging, not an official option.
    if showtrace && rem(j-1,5)==0
      fprintf('Iteration %2i: relative residual norm = %.3e\n',...
        j-1,normres(j));
    end
    
    % Give up?
    if j==maxiter+1
      warning('CHEBFUN:gmres:maxiter',...
        'Maximum number of iterations reached.')
      showtrace = false;
      break
    end
    
    % Reorthogonalize (for research only--not an official option).
    if rem(n,Inf)==0,  [Q,Rtemp] = qr(Q,0);  end
      
  end   % inner iterations
  
  y = R(1:n,1:n)\(P(:,1:n)'*QTb);                % least squares soln
  u0 = Q(:,1:n)*y;                               % new part of solution
  if ~isempty(M2inv), u0=M2inv(u0); end
  u = u + u0;                                    % solution
  v = L(u0);
  if ~isempty(M1inv), v = M1inv(v); end
  r = r - v;                                     % new residual
  if showtrace, fprintf('  (restart)\n'), end
end   % outer iterations

if j >= maxiter
    flag = 1;
end
if nargout < 2
    fprintf('\n  Final relative residual: %.3e\n\n',normres(end))
end

end   % main function

